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Abstract 

This paper is devoted to the numerical approximation of a degen- 
erate anisotropic elliptic problem. The numerical method is designed 
for arbitrary space-dependent anisotropy directions and does not re- 
quire any specially adapted coordinate system. It is also designed to 
be equally accurate in the strongly and the mildly anisotropic cases. 
The method is applied to the Euler-Lorentz system, in the drift-fluid 
limit. This system provides a model for magnetized plasmas. 
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1 Introduction 

This paper discusses the numerical resolution of degenerate anisotropic ellip- 
tic problems of the form: 



where Ocl 2 orR 3 , f £ is a given function, b is a normalized vector field defin- 
ing the anisotropy direction and e measures the strength of this anisotropy. 
In this expression V and V- are respectively the gradient and divergence op- 
erators. The unit outward normal at x G dQ is denoted by v. In the context 
of plasmas, e is related to the gyro period (i.e. the period of the gyration 
motion of the particles about the magnetic field lines), and the anisotropy 
direction b satisfies b = B/\B\ with the magnetic field B verifying V ■ B = 0. 
Eq. (11. ip may also arise in other contexts, such as rapidly rotating flows, 
shell theory and may also be found when special types of semi-implicit time 
discretization of diffusion equations are used. 

The elliptic equation is not in the usual divergence form due to an ex- 
change between the gradient and divergence operators. However, the method- 
ology would apply equally well to the operator V • ((6(8) b) ■ V0)), up to some 



(&• V) (V-(60 e )) + e0 e 
{b- v) V- {b<jf) 



f, in n, 

on <9f2, 



(1.1) 
(1.2) 
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simple changes. The expression considered here is motivated by the appli- 
cation to the Euler-Lorentz system of plasmas. This application has already 
been considered in a previous study [13] but we introduce two important 
developments. First the present numerical method does not request the de- 
velopment of a special coordinate system adapted to b. In [13], b was assumed 
aligned with one coordinate direction. Second, the present paper considers 
Neumann boundary conditions instead of Dirichlet ones as in [TJ] . Although 
seemingly innocuous, this change brings in a considerable difficulty, linked 
with the degeneracy of the limit problem, as explained below. 

A classical discretization of problem ( ll.ip . ( II. 2p leads to an ill-conditioned 
linear system as e — > 0. Indeed setting formally e = in (II. ip . (jl.2p . we get: 

-(&• V)V- {biff) = / (0) , in ft, (1.3) 
(6-i/)V-(6^) = 0, on dn, (1.4) 

with = lim e _ i ,o f e - The homogeneous system associated to (11.31) . (11.4j) 
admits an infinite number of solutions, namely all functions x/i satisfying 
V-(&0) = 0. This degeneracy results from the Neumann boundary conditions 
( II .4p and would also occur if periodic boundary conditions were used. On the 
other hand, ( II. 3p is not degenerate if supplemented with Dirichlet or Robin 
conditions, which was the case considered in [13]. A standard numerical 
approximation of ( II. 3p . (I1.4p generates a matrix whose condition number 
blows up as e — » 0, leading to very time consuming and/or poorly accurate 
solution algorithms. 

To bypass these limitations, we follow the idea introduced in [12] and 
use a decomposition of the solution in its average along the 6-field lines and 
a fluctuation about this average. This decomposition ensures an accurate 
computation of the solution for all values of e. In [12], this decomposition 
approach was developed for a uniform b and a coordinate system with one 
coordinate direction aligned with b. To extend this approach to arbitrary 
anisotropy fields b, a possible way is to use an adapted curvilinear coordinate 
system with one coordinate curve tangent to b. This is the route followed by 
[I], which proposes an extension of [12] in the context of ionospheric plasma 
physics, where the anisotropy direction is known analytically (given by the 
earth dipolar magnetic field). The approach developed here is different and 
aims at a method which does not request the generation of special curvilinear 
coordinates. Indeed, in the general case, computing such coordinates can be 
complex and costly, especially for time-dependent problems where b evolves 
in time. 
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For this purpose, we solve a variational problem for each of the terms of 
the decomposition. The main difficulty lies in the discretization of the func- 
tional spaces in which each component of the solution is searched. In the 
present paper, this difficulty is solved by introducing two kinds of variational 
systems, one corresponding to a second-order elliptic problem (for the aver- 
age) and one, to a fourth order system (for the fluctuation). An alternative 
to this method is proposed in [TO] . It avoids the resolution of a fourth-order 
problem at the price of the introdution of Lagrange multipliers which lead 
to a larger system. In the present paper, we design a method which breaks 
the complexity of the problem in smaller pieces and requires less computer 
ressources. 

As an application of the method and a motivation for studying problem 
(11.31) , fll.4p , the drift-fluid limit of the isothermal Euler-Lorentz system is con- 
sidered. These equations model the evolution of a magnetized plasma. In this 
case, the anisotropy direction is that of the magnetic field and the parameter 
e is the reciprocal of the non dimensional cyclotron frequency. The drift-fluid 
limit e — > of the Euler-Lorentz system is singular because the momentum 
equation becomes degenerate. In this paper, we propose a scheme able to 
handle both the e ~ 1 and e <C 1 regimes, giving rise to consistent approxi- 
mations of both the Euler-Lorentz model and its drift-fluid limit, without any 
constraint on the space and time steps related to the possible small value of 
e. Schemes having such properties are referred to as Asymptotic-Preserving 
(AP) schemes. These schemes are particularly efficient in situations in which 
part of the simulation domain is in the asymptotic regime and part of it is 
not. Indeed, in most practical cases, the parameter e assumes a local value 
which may change from one location to the next or which may evolve with 
time. 

The usual approach for dealing with such occurences is through domain 
decomposition: the full Euler-Lorentz model is used in the region where 
e = 0(1) and the drift-fluid limit model is used where £ < 1. There are 
several drawbacks in using this approach. The first one is the choice of the 
position of the interface (or cross-talk region), which can influence the out- 
come of the simulation. If the interface evolves in time, an algorithm for 
interface motion has to be devised and some remeshing must be used to en- 
sure compatibility between the mesh and the interface, which requires heavy 
code developments and can be quite CPU time consuming. Determining 
the right coupling strategy between the two models can also be quite chal- 
lenging and the outcome of the numerical simulations may also depend on 
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this choice. Because these questions do not have straigthforward answers, 
domain decomposition strategies often lack robustness and reliability. Here, 
using the original model with an AP discretization method everywhere pre- 
vents from these artefacts and permits to use the same code everywhere for 
both regimes. 

We conclude this introductory section by some bibliographical remarks. 
In magnetized plasma simulations, many works are based on the use of curvi- 
linear coordinate systems where one of the coordinate curves is tangent to 
the magnetic field (see e.g. [33] , the gyro-kinetic and gyro- fluid developments 
[21 [191 [22j El] and the many attempts for generating specialized coordinate 
systems [H EJ HO Q2J [23 ESI ES]). The present work, together with [TO] is 
one of the very few attempts to design numerical methods free of the use of 
special coordinate systems (see also [31]). The key idea behind this method is 
the concept of Asymptotic Preserving (AP) schemes as described above. AP- 
schemes have first been introduced by S. Jin [25] in the context of diffusive 
limits of transport models. They have recently found numerous applications 
to plasma physics in relation e.g. to quasineutrality [31 El [HI [HI E5] and 
strong magnetic fields jTUl E21 EE2] as well as to fluid-mechanical problems 
such as the small Mach-number limit of compressible fluids [16J. Other ap- 
plications of AP-schemes can be found in [3 El EDJ [23 Ell EI] . Numerical 
methods for anisotropic problems have been extensively studied in the liter- 
ature using numerous techniques such as domain decomposition techniques 
[2~Tl E7J, Multigrid methods, smoothers [3D], the hp-Unite element method 
|32j . However, these methods are based on a discretization of the anisotropic 
PDE as it is written. The method presented here as well as in [131 E21 ED] re- 
lies on a totally different concept, namely viewing the anisotropy as a singular 
perturbation and using Asymptotic-Preserving techniques. 

This paper is organized as follows. In section E] the solution methodology 
for the degenerate anisotropic elliptic problem ( II. ip . (jl.2p is detailed. Sec- 
tion E] is devoted to the discretization strategy. In section H] the drift-fluid 
limit of the isothermal Euler-Lorentz system is introduced. The AP-scheme 
is derived, giving rise to the anisotropic elliptic problem (11. ip . (II. 2p . The 
numerical method for the anisotropic elliptic problem is validated in sec- 
tion El Finally a numerical application to the Euler-Lorentz system is given 
in section El 



5 



2 A decomposition method for degenerate ani- 
sotropic elliptic problems 



We first present the methodology in the simpler case of a uniform 6-field. 
The method will then be extended to an arbitrary 6-field. 

2.1 Overview of the method in the uniform 6- field case 

A two dimensional configuration is considered in this section, with the posi- 
tion variable (x, y) belonging to a square domain (x, y) G Vt = [0, 1] x [0, 1] C 
M?. The b field is assumed uniform, equal to the unit vector pointing in the y 
direction. In this case, the singular perturbation problem (ll.ip . (II. 2p reads: 

f e (x,y), in ]0,l[x]0,l[, (2.1) 

0, fory = 0oiy= 1. (2.2) 
We assume that: 

lim ( - / f £ (x,y) dy) exists and is finite, Va; G [0, 1]. (2.3) 
s^o \e J J 

This framework is similar to [12] . Here, we recall the bases of the method- 
ology. The problem is well posed for all e > but a standard discretization 
may lead to ill-conditionned matrices when e < 1. Indeed if e is formally set 
to zero, we get the following degenerate problem 

f®(x,y), in ]0,l[x]0,l[, 

(2.4) 

0, for y = or y = 1, 

assuming that f e has the following expansion f e = /(°) + ef^ + o(e). This 
system admits a solution under the compatibility condition f^(x, y) dy = 
for all x G [0, 1], which is satisfied thanks to hypothesis (12. 3p . However the 
solution is not unique. Indeed, if ip verifies (12. 4 j) then ip + ( is also a solution 
for all functions ( = ((x) which depend on the ^-coordinate only. 



e(p £ (x,y) 



dy 2 
d_ 
dy 



<b £ (x,y) 
<j) e (x,y) 



2 



■tp(x,y) 



dy' 
dy 
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On the other hand, the limit (jr- ' = lim £ _ ) . (fi £ is unique. Indeed, it is easy 



to see that the solution ip of (12.41) such that ijj(x,y) dy = for all x G [0, 1] 
is unique. Since is a particular solution of (|2.4|) . it can be written 

0(°) =<ij, + (;(x). (2.5) 

In order to determine ( , we integrate (12. ip with respect to y and get 

i x ,i 

(p £ (x,y)dy = - f £ (x,y)dy, (2.6) 
e Jo 

Taking the limit e — > in this equation and inserting (12. 5p . we get ((x) = 
Jo f^( x i y) dy, which determines (J)^ uniquely. 

Now, if a standard numerical method is applied to (12. ip . (12. 2p . the re- 
sulting matrix will be close, when e <C 1, to the singular matrix obtained 
from the discretization of (12. 4p . Therefore, its condition number will blow 
up as e — > 0, resulting in either low accuracy, or high computational cost. 
To overcome this problem, we decompose <fi e according to 

<p e = p e + q e , P £ (x)=[ £ (x,y)dy, (2.7) 

Jo 

i.e. p e is the average of (j) e along straight lines parallel to b and q e is the 
fluctuation of the solution with respect to this average. p £ and q e satisfy: 

?f(x,y) = 0, W(x,y)eQ, (2.8) 
dy 

[ q £ (x,y)dy = 0, VxG[0,l]. (2.9) 

They are orthogonal for the scalar product of L 2 , i.e. j n p e q e dx dy = 0. 
Inserting this decomposition into (12. 6p yields 

f{x) = -J f e (x,y)dy,Vxe [0,1]. (2.10) 
Moreover, p e satisfies 



o 



limj/(x) = lim - / f £ (x,y)dy= [ 

E ^0 e ^o £ J Q J Q 



f 1 \x,y)dy = ax), 
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where ( is defined by (j2.5p . Now, q £ is the solution of the following problem: 

2 



dy 



q £ (x,y) + eq £ (x,y) = f (x, y), V(x, y) G [0, 1] x]0, 1[, (2.11) 

/ q £ (x,y)dy = 0, for xe [0,1], (2.12) 
Jo 

^-q e (x,y) = 0, for y = 0ory = l, (2.13) 

where ^ 

f = / £ - / fdy = f-ef, 



is the projection of f e on the space of functions satisfying (I2.9|) . Compared 
to (EH), (Q, system (123T]) - (1213"]) involves the additional condition (|2TT2|) . 
This condition is important: it makes the system uniformly well-posed when 
e 0. Additionally, the limit system is 

V>(s,y) = / (0) , V(x,y)G[0,l]x]0,l[, (2.14) 
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/ q (0) {x,y)dy = 0, for xG [0,1], (2.15) 
Jo 

d 

—q {0) {x,y) = 0, for y = or y = 1, (2.16) 

and has a unique solution equal to ijj. Consequently, as e — > 

e = p e + q £ -> C + "0 



Therefore, the proposed decomposition leads to two uniformly well-posed 
problems when e — > 0, which allows to reconstruct the limit solution of 
the original problem. 

The numerical approximations of conditions (I2.10p or ( 12.121) is delicate if 
the mesh is not aligned with the y coordinate axis. In order to overcome this 
problem, a weak formulation is introduced. Define V = H l (0, 1), K — {v G 
V | d y v = 0}. Then, e is the solution of the variational formulation 

Find <j) e G V such that 

^^dxdy + e [ <p £ i> dx dy = [ f e ipdxdy, G V. (2.17) 

8 



Let K L be the orthogonal space to K in L 2 (0, 1). Now, the decomposition 
( 12. 7p . corresponds to the decomposition of e on K and K^. Indeed, it is 
easily checked that p e £ K and q e £ K L and they are orthogonal, as already 
noticed. Now, inserting ip £ K in (I2.17p . we get that p e is the solution of 

Find p e £ if such that 

/ V --f £ )ipdxdy = 0, W> £ AT, (2.18) 

which means that p e is the orthogonal projection of e~ 1 f E onto if. Now, 
inserting ip £ if in ( I2.17P leads to 

Find q e £ if such that 

^-^dxdy + e f q £ ipdxdy = [ (f - ep £ ) ip dx dy, W> £ K x , (2.19) 

which is the variational formulation of (I2.1ip - (l2.13p . 

The use of these variational formulations allows for the discretization of 
( 12. ip . (12 .21) on arbitrary meshes compared to the anisotropy direction. This 
is an important advantadge over the strong formulations (I2.10p or (12. 11ft - 
( I2.13P . These formulations are now generalized to arbitrary anisotropy fields 
b in the next section. 



2.2 Presentation of the method for a general anisotropy 
field 

2.2.1 Preliminaries 

This subsection is devoted to the resolution of degenerate elliptic problems 
( II. ip . ( II. 2p for general anisotropy fields b. we first introduce the space 

V = {0 £ L 2 (fi) /V • (60) £L 2 (fi)}, 

K = {0 £ V/V • (60) = onfi}, 

W = {heL 2 (Q)/(b-V)heL 2 (Q)}, 

W = {heW /(b-v)h = Oon dtt}. 

The projection of a function on K is the generalization of the average op- 
eration ( 12.1 Op . while the projection on K 1 - corresponds to computing its 
fluctuation. The space Wo is used to characterize K^. The projections on 
K and K are well-defined thanks to the: 
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Theorem 2.1. We have the following properties 

1. K is closed in L 2 (Q). 

2. Wo equiped with the norm || h ||w = || {b ■ V)h ||l 2 (^) is a Hilbert space 
and {b ■ V)Wo is a closed space of L 2 {VL). 

3. K L = (b- V)W . 

Proof. 1) Let n G V such that n — > in L 2 (Q). Then, n — > in the 
distributional sense and the operation — V ■ (60) is continuous for the 
topology of distributions. Therfore, V • (60) = 0, which shows that G V. 

2) Wo is a Hilbert space for the norm || h \\ = \\ h + || (b-V)h \\L 2 (n)- Ac- 
cording to the Poincare inequality, the norms || || and || ||w are equivalent. 
The closedness of Wo for the L 2 topology follows from 3). 

3) The inclusion (b ■ V)Wo Q K ± is obvious. We sketch the proof of the 
converse inclusion and leave the details to the reader. We make the hy- 
pothesis that all 6-field lines are either tangent to a non-zero measure set 
of dQ or intersect dfl at two points x_ and x + such that ±(b ■ v)(x±) > 0. 
The points x_ and x + are called the conjugate points of the 6-field line and 
are respectively the incoming and outgoing points of this field line to the 
domain. These assumptions can certainly be weekened at the expense of 
technical difficulties which are outside the scope of this paper. Let if; G K 1 - . 
By taking the primitive of ip along the fe-field lines, there exists G W such 
that if; — (b ■ V)0. We can additionally impose that = on <9f2_ where 
dn± = {xedn \ ±(b- u)(x) > 0}. Let 9 G K. We have 

0= f i)6dx= I \b-V)(p9dx = [ {b-u)(j)6dS(x), (2.20) 
Jq Jn Jon 

where dS(x) is the superficial measure on dQ. Since, 9 G K its values at 
conjugate points are related by a linear relation. In particular, they can be 
taken simultaneously non-zero. Then, since the values of on dVL_ vanish, 
( I2.20p implies that the values of on dVt + vanish as well. Consequently, 
(b ■ z/)0 = on dQ, which shows that G Wo- This proves the result. □ 

Therefore, we can decompose e uniquely as 

£ = p £ + <f , p £ G K, q e G iT ± , (2.21) 
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and state problem (II. II) . (jl.2p as 

- (6 ■ V) (V • (bq £ )) + e{f + q £ ) = f £ , in fi, (2.22) 

(6-z/) V- (&<f) = 0, in (2.23) 

p e G /ST and q £ eK x . (2.24) 

Next, we introduce the variational approach. We multiply f)2.22p by a test 
function ijj G V, and integrate it on f2. Using a Green formula together with 
the boundary condition (12.231) . we find that 



V ■ (bq £ ) V • (&V0 dx + e / (p e + q £ )^dx = / / e ^- (2-25) 
n Jn Jn 

The aim now is to decompose problem ( 12.251) into a problem for p £ and a 
problem for q £ . Hence in the following two subsections the test function ip is 
chosen successively in K and in K^. 

2.2.2 Equation for p £ G K 

Chosing ip = r G K in fl 2 . 2 5 j) . we obtain the problem 

Find p £ G K such that / (ep £ - f £ ) r dx = 0, Vr G K (2.26) 

This problem admits a solution in K which is uniformly bounded in L 2 (Q) 
as e — > under the compatibility condition 

lim ( - / f £ rdx) exists and is finite, MrEK. (2.27) 



Assuming that f £ has the following decomposition 

r = f (0) + + o(e) . 

in L 2 (f2), this condition implies that G K^. Next, since ep £ — f £ G -fT -1 , 
according to Theorem 12.11 there exists g £ G Wo such that 

ep £ -f £ = (b-V)g £ . (2.28) 

Taking the product with b and the divergence of the result, we obtain the 
following 
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Proposition 2.1. p £ is given by 

p £ = ^-(f £ + b- Vg £ ) in a (2.29) 

where g e satisfies the problem: 

-V -((b®b)Vg £ ) = V-(fb) into, (2.30) 

(b-u)g £ = ondVl, (2.31) 

or, in variational form 

Find g £ G Wo such that 

[ {b ■ Vg £ ){b • V0) = / f £ b ■ V0 dx, V0 G W . (2.32) 
Jn Jn 

2.2.3 Equation for q £ G K L 

Taking %p = s G K 1 - in ( I2.25|) gives: 

/ V • (bq £ ) V • (6s) dx + e f q £ sdx= [ f £ sdx. (2.33) 
Jn Jn Jn 

But since q £ and s G K 1 - , theorem 12 . 1 1 implies that there exists h £ and 9 G Wo 
such that q £ = b ■ Vh £ and s = b ■ V0. Therefore, we get the following 

Proposition 2.2. q £ is given by: 

q e = b- Vh £ , (2.34) 
where h £ satisfies the following fourth-order problem: 

- V ■ [(& ® b) V(V ■ (b <g> b)Vh £ )} + eV ■ ((b ® &)V/i £ ) = V ■ (bf £ ), in Q, (2.35) 

(6 ■ i/)V • ((6 ® 6)V/i e ) = 0, on (2.36) 

(&. z,)/i £ = 0, ondVL, (2.37) 

or, m variational form 

Find h £ G Wo such that 

[ V • ((& ® &)V/i £ ) V • ((& ® &)V0) dx + e [ (b- Vh £ ) (b ■ V0) = 

/ £ (6 • V0) rfx, (2.38) 
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The resolution of problem (11. ip . ( II. 2p can be summarized in the following 
proposition. 



Proposition 2.3. If f £ satisfies (2.21), problem M.l\) . M.2\) is formally 



equivalent to the two problems h2.29\) . h2. 3U\) . \2.31\) on the one hand and 

( CT , mm, wm, ( fOTp . 

Remark 1. In ITUf . the characterization of K as (b ■ V)Wo is not used. 
Instead, the constraint that q G K is taken into account through a mixed 
formulation. The number of unknowns and the size of the problem are there- 
fore larger in fJTj]/ than in the present work. In practice, the resolution of 
the fourth order problem \2. 35\) . 112.36}) . \2.31ty can be reduced by solving two 



second-order problem, as shown below. Therefore, the introduction of a fourth 
order problem does not bring specific difficulties. 



2.2.4 Extension to non-homogeneous Neumann boundary condi- 
tions 

The application targeted in this paper, and detailed in section HJ requires 
the handling of non-homogeneous Neumann boundary conditions. In this 
subsection <j) £ is solution to the following inhomogeneous Neuman problem: 

e<j) £ - (b- V)(V- (b(j) £ )) = b-VK + fl, on ft, (2.39) 
(b- v)V ■ (b(j) £ ) = -(b-v)K, on dVl. (2.40) 

where n is a given function in W. We denote by f\ — b ■ Vk and by f £ = 

fi + /!■ 

Using the same decomposition (I2.2ip as before, we find that p £ satisfies 
QXHD and g £ is the solution of (1230]) . (I2T3TD or with f £ replaced by 

/| (and satisfying (12.271) ). Similarly, q £ satisfies (12.341) and h £ is the solution 
of (I2~35|) . (I236j) . ( I237j) . or of (l2T38|) with '0' at the right-hand side of ( l2~36|) 
replaced by (b ■ u)k, the other terms being unchanged. The details are left 
to the reader. 



3 Space discretization 

The problem is discretized using a finite volume method. The domain is 
decomposed into a familly IZ of rectangles Mj_ 1 / 2 j-i/2 =]xi-i, ^i[ x ]%'-i' Vjl 
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with Xi = iAx and yj = jAy. We look for a piecewise constant approximation 
p\ of p e on each Mj_i/ 2 ,j-i/ 2 and denote by Pi-\/2,j-x/2 its constant value on 
this rectangle. The function g e is approximated by a constant function on 
a dual mesh V, consisting of rectangles T> it j =]xj_i/ 2 , £i+i/2[ x ]yj-i/2, Uj+i/2[ 
where Xj_i/2 = (i — 1/2) Ax, yi-1/2 = (i — l/2)Ay. Then g £ is approximated 
by a piecewise constant function <?|> with its constant values denoted by gfj. 
We approximate (I2.29P by 

= - (rC^i,^.!) + 6(a; i _i ) y i _i) • W) Hj _ 4 ) . 

We now define approximations (6 ■ V) app and V ■ ( ■ 6) apP of operators 
^ t— ?• (6 • V^) and $ t— > V ■ (b $) such that they are discrete dual operators 
to each other. For this purpose, we define L-ji and Ld the space of piecewise 
constant functions on meshes of types 1Z and T> respectively. 

Definition 3.1. The operator (b ■ V) : L-jy — > L-ji is defined by 

= ■ ( (, 2Ax + 2A^ J' (3.1) 



2Ay 2Ay 
T/ie operator V • ( • 6) app -' — ► £x> defined by 

(V ■ (&$)app)i,j = 



2 Ax 

1 1 \ 



2 ' j 2 



2Aa: v '-a'"'-*' 2 Ay 

L ; ft.(x H) y i+ i) - —Ux^y j+k ) ) ti 3+ 



2Ax v 4 "2'^+2' 2Ay 



(3.2) 



Proposition 3.4. (b ■ V) and V • ( b ■ ) are adjoint operators to each 

1 \ / app V / app J r 
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Proof. Easy and left to the reader, thanks to a discrete Green formula. □ 



Next, we define (V ■ ((b <g> b) ■ V)) by the composition of the two operators 
(b ■ V) and V • (b ■ ) : 

\ / app V / app 

Definition 3.2. We define: 

(V • (b ® b ■ W)) app = (V • (• b)) app o (b ■ W) app , (3.3) 
where o is the composition operation. 

Finally, the approximation of problem (I2.30p . (I2.3ip is by solving the discrete 
problem for the piecewise constant function g on T>: 

(V-(6(®)6-V)) app ^ = (V-(6/)) app , (3.4) 

together with Dirichlet boundary conditions on g, where / is a piecewise 
constant function on 1Z. 

Now, problem (12.351) . (I2.36p . (I2.37P for q £ can be decomposed in two 
decoupled second-order elliptic problems of the type (12.301) . (I2.3ip and can 
be solved by a similar method. Indeed by setting u = —V • ((b (g) 6)V/i), 
we get that (12.351) . (12.361) . ( I2.37P is equivalent to the following two elliptic 
problems: 

V ■ ((b® b)Vu) -eu = V ■ (bf) in tt (3.5) 
(b-u)u = on dn (3.6) 

and 

- V • ((b (8) b) Vh) = u in Vt (3.7) 
(b-v)h = on dtt. (3.8) 

To summarize, the resolution of problem (11.11) . (11.21) reduces to three 
independent resolutions of problems similar to (13.41) . 
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4 Application to the Euler-Lorentz system in 
the drift limit 



4.1 Introduction 

In this section the drift-fluid limit of the isothermal Euler-Lorentz is investi- 
gated. This regime is representative of strongly magnetized plasma, for which 
the pressure term equilibrates the Lorentz force. It is obtained by letting a 
dimensionless parameter e, representing the non-dimensional gyro-period as 
well as the square Mach number, go to zero. This limit is singular because 
the momentum equation in the direction of the magnetic field degenerates. 
Since the field may not be uniformly large, we wish to derive an Asymptotic- 
Preserving (AP) scheme which guarantees accurate discretizations of both 
the limit regime for strongly magnetized plasma (e < 1) and the standard 
Euler-Lorentz system when the field strength is mild (e ~ 1). With this aim, 
the Euler-Lorentz system is discretized in time by a semi-implicit scheme. 

This scheme has already been studied in [13] for a uniform and constant 
magnetic field aligned with one coordinate and for physically less meaning- 
ful Dirichlet boundary conditions. The present methodology allows us to 
investigate the case of non-uniform magnetic fields and Neumann boundary 
conditions. Indeed, the anisotropic elliptic equation (II. ip . (II. 2p appears as 
the central building block of the scheme, which allows for the computation of 
the field- aligned momentum component. In this presentation, we will mainly 
focus on this aspect, the other ones being unchanged compared to 



4.2 The Euler-Lorentz model and the AP scheme 

The scaled isothermal Euler-Lorentz model takes the form: 

d t n £ + V-(n £ u £ ) = 0, (4.1) 



d t (n £ u £ ) + V • (n £ u £ <g> u £ 



TVn £ = n £ {E + u £ x B) , (4.2) 



where n £ , u £ and T are the density, the velocity and the temperature of 
the ions, respectively. Here, the electric field E and the magnetic field B 
are assumed to be given functions. The parameter e is related to the gyro- 
period of the particles about the magnetic field lines, and simultaneously to 
the squared Mach number. We refer to [13] for more details on the model, 
the scaling and the drift-fluid limit e — > 0. 
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Now we introduce the time dicretization of the model. Let B m be the 
magnetic field at time t m , \B\ m its magnitude and b m = B m /\B\ m its direc- 
tion. For a given vector field v , denote by and (v)™ its parallel and 
perpendicular components with respect to b m ie 

v = ( v )jp b m + (v)™, {v)f = v ■ b m , (tOT = b m x (v x & m ). 

Similarly, we denote by V™ and V™- the parallel gradient and divergence 
operators respective to this field. The time semi-discrete scheme proposed in 
[T5] is as follows: 

Definition 4.3. The AP scheme is the scheme defined by: 



n m+l _ n m 



At 



V ■ (nu) m+1 = 



{nu} 



\ m+l 



rax 



At 



+ V • (rax <S> u) 



T (Vn # ) 



m+l 



= n -E m+1 + (mx) m+1 x £T +1 , 
where (Vn*) m+1 groen fry, 

(Vn # ) m+1 = (Vn m )™ +1 + (Vn m+1 )7 +1 b m+1 



(4.3) 



(4.4) 



(4.5) 



By considering the scalar product of (14. 4p with b m+1 , we get 



/ \m+l / \m 

mxx) — (nu) , 1 

ri - — r-^ — — + V ■ (n e n £ (8) u e ) m ) ■ 6 m+ 



At 



= _TV rn+1 n m ■ b m+1 + n m E m+1 ■ b m+1 

and after easy computations [13], we find that (mx)™ +1 satisfies the following 
anisotropic elliptic problem: 



At 



rax 



-T At V™ +1 ( V™ +1 • ((mx) m+1 ) 



] \ m+l 

II 



TAtV^ +1 (v ■ ((rax) m+1 ) 



_l_l^m+l 



TV^ +1 n m 



+ 



LAt 



nw) m — e ( V • (rxxx ® u) 



II 



m jrim+1 



(4.6) 



n m E 



m+l 
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By setting 

1 



and by taking / — fx + / 2 with 



h 

h 



At 



b- V(V- (nu™ +1 )), 



:{nu) m - 



IT (At) 



— (b-Vri 
AV 



T At 



V • (nu®u) rn + n m E 



m rpm+1 



m+1 



(4.7) 
(4.8) 

(4.9) 



this problem can be put in the framework of (11.11) . In [13], because b was 
chosen parallel to one of the coordinate direct discretization of H4. 61) 

using finite differences could be performed. Here, for an arbitrary anisotropy 
direction b, we use the method developed in the previous sections. We do 
not detail the description of the discretization of the other equations, since 
it follows [TBI. 



The right-hand side (I4.9[) can be decomposed as /| = ffl +sf^ with 
corresponding to the first two terms and f%, to the last two one. Moreover 
if we suppose that 



n m E m+l 



m+1 1 

- — (&■ Vn m ) e K~ 

II AV ' 



(4.10) 



the compatibility condition (I2.27P is satisfied. This property amounts to 
saying that the integrated force along a magnetic field line is zero. If the 
property is not satisfied, parallel velocities of order 0(1 je) are generated, 
which is physically unrealistic (because collisions will ultimately slow down 
the plasma ions). Therefore, assuming (14.101) is physically justified. 

As in [13], we will compare the AP scheme with the classical semi-discrete 
scheme for the Euler-Lorentz model, given by: 



Definition 4.4. The 'classical' semi- discrete scheme is defined by: 

,m+l 



n 



n 



At 



+ V 



[nu\ 







/ \m+l / \m 

[nu] — (nu) 



At 



V • (nu C*D u) 



+ T (Vn) 



n m E m+l + ( nu )^ x B 



m+1 



?m+l 



(4.H) 



(4.12) 



In [13], it is shown that this scheme is not uniformly stable with respect to 
e and so that it cannot be AP. 
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Except from the parallel momentum equation, which has just been dis- 
cussed, the other equations of the model are discretized following [13]. For 
the sake of brevity, we will not reproduce their presentation here. 



4.3 Boundary conditions 

The following boundary conditions are set up for test purposes only. We 
impose Dirichlet boundary conditions on the density n m+1 = ub with ng in- 
dependent of time. For the perpendicular momentum, we impose the relation 
obtained after taking the limit when e — > in (14. 2p . 



nu™ +1 = - , 1 - b x (T Vn m + n m E m+1 ) . 

Ig m+1 \ / 

By considering the mass conservation equation at the domain boundary, we 
have 



n m+l _ n r 



At 



+ V- (bn m+1 ul +1 ) + V • = 0, on dQ. 



Therefore, as the density satisfies Dirichlet boundary conditions with time- 
independent Dirichlet values, we get 

(&• u)V ■ (bnu™ +l ) = -(&• i/)V ■ (n< +1 ), on dQ. 

Therefore, nu™ +1 is a solution to the anisotropioc elliptic problem with inho- 
mogeneous Neumann boundary conditions (12.391) . (12.40(1 . with k = —(b-i')'V- 
(nu™ +1 ). Then, we can apply the framework of section 12.2.41 When nu\\ has 
been calculated, an approximation is employed in order to provide values of 
nun in a layer of fictious cells surrounding the boundary, by using homoge- 
neous Neumann boundary conditions. The values in the fictitious cells are 
then useful to compute gradient terms which occur in the other equations of 
the Euler-Lorentz model. 



5 Numerical results for the elliptic problem 
5.1 Introduction 

In this section the efficiency of the numerical method introduced in sec- 
tions [2] and [3J for the singular perturbation problem ( 11. ip . 1 11. 2ft is investigated 
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through numerical experiments. These experiments are carried out on a two 
dimensional uniform Cartesian mesh. Two sets of test cases are presented. 
In the first one, the anisotropy, or magnetic field, is oblique, which means 
that it is assumed uniform in space, but not necessarily aligned with any 
coordinate axis. In the second set, the field direction is non uniform. In 
both cases, the strength of the anisotropy is assumed uniform and is given 
by the value of e. An analytical solution <p a is constructed for the singular 
perturbation problem (11. ip . (11. 2p and is compared with its approximation 
cf) h computed on the mesh. For the test cases, the following L 1 , L 2 and L°° 
norms are used to estimate the errors between the numerical approximation 
cj) h and the analytical solution <j) a : 



Ua~< 




L 1 


^tjlM^Vj) -<i> h (iJ)\ 


Ua\ 






Eij \M%i,yj)\ 


Ua~< 


f\ 


w _ 


(Eul^fo,v;)-#il 2 )* 


Hal 








Ua~< 


A 




maxij \<p a {x i ,y j ) 



(5.1) 



Ua\\L°° maxij |0 a (x i ,?/ j )| 



5.2 Numerical results for an oblique magnetic field 
5.2.1 Introduction and test case settings 

For these numerical experiments the simulation domain is the square Q = 
[0,1] x [0,1]. The magnetic field is defined by B = (sin a, cos a, 0), with a 
the angle of the b-field with the x-axis ranging from to n/2. In order to 
validate the numerical method an analytical solution denoted <j) a for problem 
(II. ip . (II .2p is constructed. It is written 

a (x, y) = sin (x sin(a) — y cos(a)) + b ■ VH(x, y) , 
f a (x, y) = -b ■ V(V • ((6 b)VH(x, y))) 

+ e ( sin (x sin(a) — y cos(ct)) + b ■ VH(x, y)) , 

H (?, y) = {{x - l)(y - l)xy) 3 . 

The function a is the solution of problem (II. ip . (II. 2p with the right-hand side 
fl- (p a presents itself as decomposed into p e (first terms) and q £ (second term). 
Note also that fl can be decomposed as f £ a = fa +efa with = —b-Vh 
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and h = V ■ ((&®&)V H(x, y)). The function h verifies homogeneous Dirichlet 
boundary conditions on the domain boundaries, which implies, according 
to theorem 12.11 that e K 1 - and the compatibility condition (I2.27P is 
satisfied. However, for the simulations carried out below, the construction 
of the right-hand side f e a is performed using the discrete operators (6 ■ V) app 
and V • ( ■ fr) app in order to ensure that the compatibility condition (I2.27P is 
satisfied by the discrete operators, namely /i £ ^app> where 

i^app = {0 / V • (60) app = 0} , K^ pp = (b • V) app (Wo) . 

5.2.2 Homogeneous Neumann boundary conditions 

This simulation is run with a = it/ 3. On figured], we represent the relative 
errors as functions of the mesh sizes for different values of e ranging from 10 -3 
to 10~ 9 . The curves of figured] are plotted using logarithmic decimal scales. 
We observe a linear decrease of the errors with vanishing mesh sizes, with a 
slope equal to 2, which proves that the global scheme is second order accurate. 



More importantly, we observe from figures 1(a) and 1(b), that the precision 
remains the same while e is decreased by three orders of magnitude. However, 
for the more refined grids using the smallest value of e of this simulation set 



(1CT 9 , see figure 1(c)), a slight degradation of the convergence is observed 
for small mesh sizes. 

This slight degradation can be explained. Indeed, p £ is given by a stiff 
problem, since ep e is obtained as the difference of two quantities scaling as 
e° = 0(1) (see (12T29]) . (B3Djl ). To investigate the influence of e on the 
accuracy of the approximation of p 6 , the L°° norm of the relative error made 
on p e and on V • (bp £ ) as functions of e are plotted on figure [2] 



Figure 2(a) shows a linear behavior of V ■ (bp £ ) with vanishing e (in log 
scale). To explain this feature, we note that the discretization of the second 
order operator in ( I2.30P provides a computation of e (V • (■ &)) app (p e ) with 
the precision of the linear system solver used for the computation of g £ , which 
is limited by round-off errors. This error is amplified after multiplication by 
the factor 1/e. This analysis still holds for the accuracy of p e as a function of 



e represented on figure 2(b) with slight differences. For the largest values of 
s, we observe a plateau (red dashed line) explained by the discretization error 
of the discrete operators. The space discretization introduced here is second 
order accurate, i.e. is 0(h 2 ) where h = max( Ax, Ay) . Since the right-hand 
side is well prepared this error only applies to the ef^- 1 ' part of f £ and is then 
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Error w.r.t the step size in space. 



Error w.r.t the step size in space. 





-1.5 



-1.4 



Size step in log scale 

(a) s = 10~ 3 . 

Error w.r.t the step size in space. 




Size step in log scale 



(c) 



10- 



Size step in log scale 



(b) e = 10- 



Figure 1: Oblique magnetic field 
test case with a = n/3: error 
norms, defined by (15.11) . for the so- 
lution d> e clS 8b function of the mesh 
size, in decimal logathimic scales, 
and for different values of e. 



proportional to eO(h 2 ) in b- Vg e , giving rise to a 0(h 2 ) consistency error for 
p e . The value of the plateau is thus only dependent of the mesh sizes and 
does not depend on the values of e. With vanishing values of e the round-off 
errors due to the linear system solver grow linearly (in log scale) until they 
reach the consistency error (0(h 2 )). This occurs for a value of e which, for 
this test case, can be estimated as approximately e = 10 -9 . For smaller e, the 
discretization error is negligible compared to the round-off errors amplified 
by the factor 1/e and the accuracy of p £ deteriorates linearly with vanishing 
e. 

The accuracy of the approximation of p e can be made totally independent 
of e under the assumption that = 0. In this case, both b ■ Vg e and f £ 
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£ (decimal log scale). 

(a) Infinity norm for V • (bp £ ) as a function 
of e in decimal log. scales. 




-20 -16 -16 -14 -13 -10 -8 -6 -4 -2 

£ (decimal log scale). 

(b) Relative error in infinity norm for p e 
as a function of e in decimal log. scales. 



Figure 2: Oblique magnetic field test case for a = it/ 3 and Ax = Ay = 1/60. 
Approximation of the p £ part of the solution. 



scale as e, providing then an approximation of p £ independent of e. The 
numerical methods introduced in [10J [12] have been developed under this 
assumption that = 0. The present paper is developed under a weaker 
hypothesis, required by the application to the Euler-Lorentz model in the 
drift-limit. This explains why a comparable accuracy cannot be reached. 
Therefore, strictly speaking, our scheme is AP for the computation of p £ 
only when = 0, or, when /(°) ^ 0, only if the round-off errors brought 
by the linear system solver are smaller than the discretization error. Still, 
it is AP without any restriction for the computation of q e (i.e. even when 

The next simulation is aimed at investigating whether the accuracy de- 
pends on the angle between b and the coordinate axes. For this purpose, 
simulations are carried out on a mesh composed of 40 x 40 cells and for a 
ranging form to tt/2. When a = the b field is aligned with the x-axis and 
when a = n/2, it is aligned with the y-axis. The relative errors are displayed 
as functions of a on figure [3] We observe that the variations of the errors are 
small on the whole range of angles. This confirms that the method provides 
accurate results, even when the mesh is far from consistent with the 6-field 
direction. 
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Error in function of the direction of the magnetic field. 



O error in L1 norm 

error in L2 norm 

+ error in infinite norm 




1 1 1 1 1 1 1 1 1 1 

0.2 0.4 0.6 0.8 1 1 .2 1 .4 1 .6 

Direction of the magnetic field. 



Figure 3: Oblique magnetic field test case for e = 10 -9 and Ax = Ay = 1/40. 
Norms of the relative error (15. ip as a function of the angle of the magnetic 
field with the x-axis a. 

5.2.3 Inhomogeneous Neumann boundary conditions 

We remark that <p £ (x, y) = 2x 2 + y 2 is an analytical solution of system (I2.39p . 
(I2.40p for f2{x, y) = e(2x 2 +y 2 ) and n = — V-(6/). For this analytical solution 
and e = 10~ 9 , we have checked that the relative error does not exceed 10" 13 . 

5.3 Numerical results for a non uniform magnetic field 
5.3.1 Introduction and test case settings 

In this subsection Q =]1, 2[x]l, 2[ and the magnetic field is given by: 

B=\B\b, & = (frin(0),-cos(0)), tan(0) = -. (5.3) 

x 

For this case, an analytical solution of (II. ip . (II. 2ft can be found. We consider 
H var defined on [1,2] x [1,2] by H var (x,y) = (1 - x) 3 (l - yf{2 - xf{2 - yf. 
According to Theorem 12. 1[ b ■ V H var £ K L . So = 1 + b ■ V H var is the 
solution of (II. ip . Q1.2j) when the right-hand f £ of (ll.ip has the expression 

f = -b ■ V (V • (b ® 6) VH var ) +e(l + b- VH var ) . 
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5.3.2 Homogeneous Neumann boundary conditions 



On figures 4(a), 4(b) and 4(c), we have represented the relative errors as 
functions of the mesh size when e goes from 10~ 3 to 10~ 9 . We observe that 
all the three norms decrease when the mesh sizes decrease, in a similar fashion 
as in the oblique uniform 6-field. 



the step size in space. 




-1.6 -1.5 -1.4 -1.3 -1.2 

Size step in log scale 



(a) e = 1(T 



Error w.r.t the step size in space. 




Size step in log scale 



-4.4 
-4.6 



c -5.2 
£ -5.4 



Error w.r.t the step size in space. 




Size step in log scale 



(b) e = 10- 



Figure 4: Non uniform magnetic 
field test case: error norms, de- 
fined by (I5.ip . for the solution <p £ 
as a function of the mesh size, in 
decimal logathimic scales, and for 
different values of e. 



(c) £ = 10- 



5.3.3 Inhomogeneous Neumann boundary conditions 

We take the test case of subsubsection 15.2.31 again, and we find a similar 
conclusion: with e = 10~ 9 , the relative error in L°° norm does not exceed 

io- 11 . 
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6 Numerical results for the Euler-Lorentz sys- 
tem in the drift limit 



6.1 Introduction and test case settings 

This part is devoted to the validation of the AP-scheme (14. 3 j) , (14. 4p , (14. 5 p for 
the Euler-Lorentz system. Due to the lack of analytical solutions, the valida- 
tion procedure consists in comparisons of the AP-scheme with the classical 
discretization (I4.12p . The classical discretization is subject to a CFL stabil- 
ity condition that imposes the time step to resolve (i.e. to be smaller than) 
the fastest time scales involved in the system. These time-resolved simula- 
tions require a time step which scales like \fe (because the CFL condition 
involves the acoustic wave speed which scales like 1/y/e). The AP-scheme is 
designed to be stable independently of e when e — > 0. In these situations, 
the time step cannot resolve the fastest time scales involved in the system, 
which leads to under-resolved simulations. The stability of the AP-scheme 
in under-resolved situations has be demonstrated in [13]. In this case, the 
requested CFL condition only involves the fluid velocity, which is an 0(1) 
quantity, and not the acoustic speed [13] and explains the possibility of using 
large time steps, independent of e. We want to check this feature again when 
the scheme is equipped with our new elliptic solver. 

Two test cases are presented, one for an oblique uniform magnetic field, 
another one for a non uniform magnetic field with the same expressions as 
in section [51 In both cases, the electric field is chosen as E = (0, 0, B x + B y ), 
where B x and B y are the components of the magnetic field. The initial 
condition is defined by the following uniform data: n — 1, (nu) x = 1, {nu) y = 
— 1 and (nu) z = which defines a stationary solution of the Euler-Lorentz 
system. A local perturbation of order e in then applied to this stationary 
state and the evolution of the system is observed for both the AP and the 
classical schemes. 

6.2 Numerical results for an oblique uniform magnetic 
field 

The results for the AP and the classical schemes are compared on figure [5] 
in a resolved case. Both schemes provide comparable results. However we 
observe the formation of a thin boundary layer on the domain frontiers for 
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the AP-scheme but it is not responsible for the development of an instability. 

Next we consider the same test case with an under-resolved time step At 
which is 10 times larger than the time step provided by the CFL condition 
of the classical scheme. These simulation results are collected on figure [6j 
In this case, the conventional scheme leads to unstable results contrary to 
the AP scheme and proves the capability of the AP-scheme to provide stable 
computations for time steps that resolve neither the acoustic wave-speed nor 
the gyration period. 



6.3 Numerical results for a non uniform magnetic field 



For the non uniform case, n = 1, (nu) x = 1, (nu) y = —1 and (nu) z = are 
not stationary solutions to the Euler-Lorentz system. In particular, with the 
chosen initial condition, sharp boundary layers are generated. But the the 
AP scheme can still be compared with the classical scheme in the resolved 
case for a validation procedure. Then we take the same initial conditions as 



7(f) , 7(e)) 



for the oblique magnetic field case. Figures (7(b) , 7(a) , 7(d) , 7(c 
show that the two schemes provide similar results. 

Next we consider the under-resolved time step lOAi. In this situation 
show that the classical scheme is unstable. By contrast, 
demonstrate that the AP-scheme provides stable results. 



Fig. 


8(a) 


8(c), 


8(e) 


Fig. 


8(b) 


8(d), 


8(f) 



The increased numerical diffusion generated by the large time step gives rise 
to a widening of the boundary layer. Keeping the boundary layer accurate 
would require some mesh refinment in the vicinity of the boundary. This 
point is deferred to future work. 

Moreover as the initial conditions of the present test case are not sta- 
tionary solutions of the Euler-Lorentz model, it is important to check if the 
results obtained in the non resolved case by the AP scheme correspond to the 
proper limit regime. So we compare the results obtained with and without 
the local perturbation on the initial conditions. The difference between the 
results obtained with the two simulations remain of the same order as the 



perturbation of the initial condition. Fig. (9(a) 9(b) , 9(c) ) present the differ- 



ence between the solutions obtained with the perturbed and non-perturbed 



initial condition, for n, (nu) 



\nu) 



The figures show that this difference is 



actually of 10 -10 for the density and 10 -6 for the momenta. The difference 
with the value of e — 10 -9 , can be explained by the accumulation of the 
truncation error over the simulation time. 
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7 Conclusion and perspectives 



A numerical method for degenerate anisotropic elliptic problems has been 
investigated. This method is based on a variational formulation together 
with a decomposition of the solution. This problem has been applied to 
the resolution of an Asymptotic-Preserving scheme for the isothermal Euler- 
Lorentz system. Numerical simulations demonstrate the ability of the scheme 
to handle under-resolved situations where the time-step exceeds the CFL 
stability condition of the classical scheme. 

Forthcoming works will be devoted to the generalization of this approach 
for the full Euler system with a non linear pressure law. In this case non 
linear anisotropic elliptic problem have to be handled. Moreover we can also 
deal with the more physical situation of a plasma constituted by a mixture 
of ions and electrons. In this situation the model can be described by the 
two-fluid Euler-Lorentz system coupled with quasi-neutrality equation. 
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(a) n (AP-scheme) . 



(b) n (classical scheme). 





(e) nu y (AP-scheme). (f) nu y (classical scheme). 

Figure 5: Euler-Lorentz test case for an oblique magnetic field in the resolved 
case at time t = 3.95 10 -6 s: density (n) and momentum (nu x , nu y ) com- 
puted by the AP-scheme (left) and the classical scheme (right) for e = 10 -9 
and Ax = Ay = 1/40. The angle of the magnetic field with the x-axis is 
a = 7r/3. 
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(a) n (AP-scheme). 



(b) n (classical-scheme). 





(e) nu y (AP-scheme) (f) nu y (classical-scheme) 

Figure 6: Euler-Lorentz test case for an oblique magnetic field in the under- 
resolved case at time 3.95 10~ 5 s: density (n) and momentum (nu x , nu y ) 
computed by the AP-scheme (left) and the classical scheme (right) for e = 
10 -9 and Ax = Ay = 1/40. The angle of the magnetic field with the x-axis 
is a = 7r/3. 
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(a) n (classical-scheme) 



(b) n (AP-scheme) 




(c) nu x (classical-scheme) (d) nu x (AP-scheme) 




(e) nu y (classical-scheme) (f) nu y (AP-scheme) 

Figure 7: Euler-Lorentz test case for a non uniform magnetic field in the 
resolved case at time t = 3.95 10~ 6 s: density (n) and momentum (nu x , 
nu y ) computed by the AP-scheme (right) and the classical scheme (left) for 
s = 10~ 9 and Ax = Ay = 1/40. 
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(a) n (classical-scheme) 



(b) n (AP-scheme) 





(c) nu x (classical-scheme) 



(d) nu x (AP-scheme) 





(e) nu y (classical-scheme) 



(f) nu y (AP-schcmc) 



Figure 8: Euler-Lorentz test case for a non uniform magnetic field in the 
under-resolved case at time t = 3.95 10 _5 s: density (n) and momentum(n« I , 
nu y ) computed by the AP-scheme (right) and the classical scheme (left) for 
e = 10~ 9 and Ax = Ay = 1/40. 
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(a) Difference for n 



(b) Difference for nu. 



x^o B 




Figure 9: Difference between the solu- 
tions obtained with an initial pertur- 
bation of order e = 10 -9 and the solu- 
tion without any perturbation for the 
variable magnetic field after 1.58 s of 
simulation in the non resolved case. 



(c) Difference for nu 
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rho with the classical scheme 
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